\subsection{Cell Transmission Model and its Properties}\label{subsec_ctm}
Our traffic model of choice is the Cell Transmission Model (CTM)
\cite{daganzo94, daganzo95}.
CTM is popular for its flexible use in macroscopic simulation.
Compared with microscopic simulators, it requires negligible
computational effort.
It can be extended to road networks (\cite{strada96a})
and urban roads with signalized intersections (\cite{lo01, almasri05}).
It appears well-suited for calibration with point measurements of aggregate
traffic variables that are routinely available (\cite{lin95, munoz04}).
It is straightforward to formulate questions of ramp metering
(\cite{daganzo93,zhang96,gomes06}) and dynamic traffic assignment
(\cite{strada96b,ziliaskopoulos00}) by posing optimization problems within CTM.

\begin{figure}[htbp]
\centerline{
\includegraphics[width=5in]{figures/fwcells.eps}}
\caption{Freeway divided into $N$ cells.}
\label{figfwcells}
\end{figure}

The freeway is divided into $N$ cells, each with one on- and one off-ramp,
indexed $1, \cdots , N$ as shown in Figure \ref{figfwcells}.
Cell $i$ is upstream of cell $i+1$.
There are two boundary conditions.
Free flow prevails downstream of cell $N$; upstream  of the freeway is an
``on-ramp'' with an inflow of $r_0$.
The flow accepted by cells $1$ is $f_0(t)$ vehicles per period or
time step $k$.
Cell $i$ is characterized by a single state variable, its density $\rho_i$, 
so the state of the freeway is the $N$-dimensional vector
$\rho = (\rho_1, \cdots , \rho_{N})$.
Vehicle movement in a cell is governed by a triangular fundamental
diagram defined by capacity, free flow speed and congestion wave speed.
(Figure \ref{fig-fdiagram}).
\begin{figure}[ht]
\centering
\includegraphics[width=3in]{figures/fd.eps}
\caption{The fundamental diagram in cell $i$ is characterized by the capacity
$F_i$, freflow speed $v_i$ and congestion wave speed $w_i$.}
\label{fig-fdiagram}
\end{figure}
Density value that corresponds to the capacity is the critical density,
and density at which the traffic speed is $0$ is the jam density.
If the density is below critical, vehicles move at free flow speed;
if it is above critical, the cell is congested, speed is lower,
and flow from the immediately upstream cell is constrained.
Thus, the state of a freeway obeys a $N$-dimensional nonlinear difference
equation.
\begin{equation}
\dot{\rho}_i(t) =
\frac{1}{l_i}\left(f_{i-1}(t, \rho(t)) -
\bar{\beta}_i^{-1}f_{i}(t, \rho(t)) + r_i(t)\right),
\;\;\; i = 1, \cdots, N,
\label{ctctm0}
\end{equation}
with
\begin{eqnarray}
f_0(t, \rho(t)) & = & \min\left\{r_0(t),\;
w_1(\bar{\rho}_1 - \rho_1(t)),\; F_1\right\},
\label{ctctm1}\\
f_i(t, \rho(t)) & = & \min\left\{v_i\rho_i(t),\;
w_{i+1}(\bar{\rho}_{i+1} - \rho_{i+1}(t)),\; F_{i+1}\right\},
\;\;\; i = 1, \cdots, N-1,
\label{ctctm2}\\
f_N(t, \rho(t)) & = & \min\left\{v_N\rho_N(t),\; F_N\right\}.
\label{ctctm3}
\end{eqnarray}
The notation is as follows:
\begin{description}
\item[$N$] is the number of cells;
\item[$i$] is the cell index, $i = 1, \cdots, N$;
\item[$l_i$] is the length of cell $i$ in \emph{miles};
\item[$t$] is the time;
\item[$\rho_i^c$] is the critical density for cell $i$,
in \emph{vehicles per mile};
\item[$\bar{\rho}_i$] is the jam density for cell $i$,
in \emph{vehicles per mile};
\item[$F_i$] is the capacity of cell $i$,
in \emph{vehicles per hour};
\item[$v_i=\frac{F_i}{\rho_i^c}$] is the free flow speed,
in \emph{miles per hour};
\item[$w_i=\frac{F_i}{\bar{\rho}_i - \rho_i^c}$] is the congestion wave speed,
in \emph{miles per hour};
\item[$\rho_i(t)$] is the density in cell $i$, in \emph{vehicles per mile};
\item[$\rho(t)$] is the $N$-dimensional vector of densities;
\item[$f_0(t, \rho(t))$] is the flow entering the first cell
in \emph{vehicles per hour};
\item[$f_i(t, \rho(t))$] is the flow leaving cell $i$,
in \emph{vehicles per hour};
\item[$r_i(t)$] is the on-ramp flow for cell $i$,
in \emph{vehicles per hour};
\item[$\beta_i$] is the off-ramp split ratio for cell $i$,
$0\leq\beta_i\leq 1$, defines which portion of the cell outflow leaves
through the off-ramp;
\item[$\bar{\beta}_i=1-\beta_i$] is the complementary split ratio for cell $i$.
\end{description}

Split ratios $\beta_1 \cdots, \beta_{N}$ are fixed.
Assume stationary demands $r_i(t) \equiv r_i$.
Each on-ramp \emph{demand} vector $r = (r_0,\cdots , r_N)$
induces a unique equilibrium \emph{flow} vector $f(r) = (f_0, \cdots, f_N)$
calculated as
\begin{eqnarray}
f_0 & = & r_0, \label{2d} \\
f_i & = & \bar{\beta}_i\: ( f_{i-1} + r_i ), \quad 1 \le i \le N. \label{2a}
\end{eqnarray}
The function $r \mapsto f(r)$ defined  by (\ref{2d}), (\ref{2a}) is bijective.
We say that demand $r$ is \emph{feasible}
if $0 \le f_i \le F_i, \, 0 \le i \le N$;
it is \emph{strictly feasible} if $0 \le f_i < F_i, \, 0 \le i \le N$;
and it is infeasible if $f_i > F_i$\footnote{Note that equilibrium flows
$f_i$, $i=0..N$, depend only on demand and split ratios,
and hence are constant, as opposed
to flows $f_i(t)$ coming from (\ref{ctctm1})-(\ref{ctctm3}).}.

Following \cite{gomes08}, we say that $i$ is a \emph{bottleneck} cell for
demand $r$ (or induced flow $f$) if $f_i \geq F_i$.
It is shown in \cite{gomes08} that if there are $K \ge 0$ bottlenecks at 
$1 \le I_1 < I_2 \cdots < I_K \le N$,
the freeway can be partitioned into $1+K$ segments $S^1, \cdots , S^{K+1}$
comprising contiguous cells as follows:
\begin{equation}
S^1=\{1, \cdots , I_1 \}, \;
\cdots , \;
S^K = \{I_{K-1}+1, \cdots, I_K\}, \;
S^{K+1} = \{I_K+1, \cdots, N\}, \label{segments}
\end{equation}
and each segment $S^j$, $1\leq j \leq K$, can be considered separately.

Moreover, system (\ref{ctctm0})-(\ref{ctctm3}) is strictly monotone, meaning
that for any initial conditions $\tilde{\rho}(0)$ and $\hat{\rho}(0)$,
such that $\tilde{\rho}(0)<\hat{\rho}(0)$, it is true that
$\tilde{\rho}(t)<\hat{\rho}(t)$ for any $k\geq0$\footnote{For two vectors
$x, y$ in ${\bf R}^N$, write
\begin{eqnarray*}
x \le y & \Leftrightarrow& x_i \le y_i, \;\;\; i=1..N,\\
x < y &\Leftrightarrow& \; x \le y, x \neq y.
\end{eqnarray*}. }

Our goal is to estimate the state $\rho(t)$ for time $t$,
$0\leq t\leq T$ --- that is, for the time period of constant demand.
Without loss of generality, we shall consider the freeway with a single
bottleneck in cell $N$ --- in other words, the we assume that the whole
freeway consists of just one segment $S^1$.

Measurement data comes from point and/ or mobile sensors. 
Point sensors measure flows $f_i(t)$, $i=0..N$.
Mobile sensors measure traffic velocities $V_i(t)$, $i=1..N$, inside the cells.
For cell $i$, density $\rho_i(t)$, flow $f_i(t)$ and velocity $V_i(t)$
are related as follows:
\begin{equation}
V_i(t) = \frac{f_i(t)}{\bar{\beta}_i\rho_i(t)}. \label{density-flow-speed}
\end{equation}
The observation equation is
\begin{equation}
y(t) = h(t, \rho(t)),\label{ctctmoutput0}
\end{equation}
where measurement $y(t)$ is either $N+1$-dimensional vector flows $y_f(t)$
in the case of point sensors,
or $N$-dimensional vector of velocities $y_V(t)$ in the case of mobile sensors.
More specifically,
\begin{equation}
y_f(t) = \left(\begin{array}{c}
f_0(t, \rho(t))\\
\vdots \\
f_N(t, \rho(t))\end{array}\right), \label{ctctmoutputf}
\end{equation}
and
\begin{equation}
y_V(t) = \left(\begin{array}{c}
f_1(t, \rho(t)) / \bar{\beta}_1\rho_1(t)\\
\vdots \\
f_N(t, \rho(t) / \bar{\beta}_N\rho_N(t))\end{array}\right). \label{ctctmoutputV}
\end{equation}
\begin{prop}\label{prop_observability}
System (\ref{ctctm0})-(\ref{ctctm3}) is completely observable
\begin{enumerate} 
\item under flow measurement $y_f(t)$ iff all cells are in free flow:
\[ f_0(t, \rho(t)) = r_0 \; \mbox{ and } \;
f_i(t,\rho(t)) = v_i\rho_i(t), \; i=1..N, \]
or all cells are congested and the mainline demand $r_0$ cannot be satisfied:
\[ f_i(t, \rho(t)) = w_{i+1}(\bar{\rho}_{i+1} - \rho_{i+1}(t)), \; i=0..N-1, \;
\mbox{ and } \; f_N(t, \rho(t)) = F_N, \]
or upstream cells are congested, the mainline demand $r_0$ cannot be satisfied
and free flow cells can only be downstream of the congested cells:
\[ f_i(t, \rho(t)) = w_{i+1}(\bar{\rho}_{i+1} - \rho_{i+1}(t)), \; i=0..M<N, \;
\mbox{ and } \; f_i(t, \rho(t)) = v_i\rho_i(t), \; i=M+1..N, \]

\item under velocity measurement $y_V(t)$ iff all cells are congested:
\[ f_i(t, \rho(t)) = w_{i+1}(\bar{\rho}_{i+1} - \rho_{i+1}(t)), \; i=1..N-1, \;
\mbox{ and } f_N(t, \rho(t)) = F_N. \]
\end{enumerate} 
\end{prop}
To get a better understanding of this proposition,
let us study simple 2-cell example.










\subsection{Example: 2-Cell Freeway Segment}
Consider freeway segment shown in Figure \ref{fig2cells}
with the following parameters:
\begin{description}
\item[$l_1 = l_2 = 1$ mile;]
\item[$\rho_1^c = \rho_2^c = \rho^c = 100$ vpm;]
\item[$\bar{\rho}_1 = \bar{\rho}_2 = \bar{\rho} = 400$ vpm;]
\item[$F_1 = F_2 = F = 6000$ vph;]
\item[$v_1 = v_2 = v = \frac{F}{\rho^c} = 60$ mph;]
\item[$w_1 = w_2 = w = \frac{F}{\bar{\rho}-\rho^c} = 20$ mph;]
\item[$r_0 = 4800$ vph;]
\item[$r_2 = 1200$ vph;]
\end{description}

\begin{figure}[htbp]
\centerline{
\includegraphics[width=2.5in]{figures/2cells.eps}}
\caption{2-cell example.}
\label{fig2cells}
\end{figure}

Piecewise affine system (\ref{ctctmstate}) for $N=2$ has
six affine modes listed in Table \ref{tab_2cells}.

\begin{table}[ht]\begin{center}
\begin{tabular}{l|l|l}
Mode & Conditions & System Equations \\

\hline
\emph{FF} &
$\begin{array}{lll}
\rho_1(t) & \leq & \rho_1^c\\
\rho_2(t) & \leq & \rho_2^c\end{array}$ &
$\begin{array}{l}
\dot{\rho}(t) = \left(\begin{array}{cc}
-v_1 & 0 \\
v_1 & -v_2 \end{array}\right) \rho(t) + \left(\begin{array}{c}
r_0 \\
r_2 \end{array}\right)\\
y_f(t) = \left(\begin{array}{cc}
v_1 & 0\\
0 & v_2\end{array}\right)\rho(t) \\
y_V(t) = \left(\begin{array}{c}
v_1\\
v_2\end{array}\right)\end{array}$\\

\hline
\emph{FC} &
$\begin{array}{lll}
v_1\rho_1(t) & \leq & w_2(\bar{\rho}_2 - \rho_2(t))\\
\rho_2(t) & > & \rho_2^c\end{array}$ &
$\begin{array}{l}
\dot{\rho}(t) = \left(\begin{array}{cc}
-v_1 & 0 \\
v_1 & 0 \end{array}\right) \rho(t) + \left(\begin{array}{c}
r_0 \\
r_2 - F_2 \end{array}\right) \\
y_f(t) = \left(\begin{array}{cc}
v_1 & 0\\
0 & 0\end{array}\right)\rho(t) + \left(\begin{array}{c}
0\\
F_2\end{array}\right) \\
y_V(t) = \left(\begin{array}{c}
v_1\\
F_2/\rho_2(t)\end{array}\right)\end{array}$\\

\hline
\emph{CF1} &
$\begin{array}{lll}
\rho_1(t) & > & \rho_1^c\\
w_1(\bar{\rho}_1 - \rho_1(t)) & \geq & r_0\\
\rho_2(t) & \leq & \rho_2^c\end{array}$ &
$\begin{array}{l}
\dot{\rho}(t) = \left(\begin{array}{cc}
0 & 0 \\
0 & v_2\end{array}\right) \rho(t) + \left(\begin{array}{c}
r_0 - F_2\\
r_2\end{array}\right) \\
y_f(t) = \left(\begin{array}{cc}
0 & 0 \\
0 & v_2\end{array}\right) \rho(t) + \left(\begin{array}{c}
F_2 \\
0\end{array}\right) \\
y_V(t) = \left(\begin{array}{c}
F_2 / \rho_1(t)\\
v_2\end{array}\right)\end{array}$\\

\hline
\emph{CF2} &
$\begin{array}{lll}
w_1(\bar{\rho}_1 - \rho_1(t)) & < & r_0\\
\rho_2(t) & \leq & \rho_2^c\end{array}$ &
$\begin{array}{l}
\dot{\rho}(t) = \left(\begin{array}{cc}
-w_1 & 0 \\
0 & v_2\end{array}\right) \rho(t) + \left(\begin{array}{c}
w_1\bar{\rho}_1 - F_2\\
r_2\end{array}\right) \\
y_f(t) = \left(\begin{array}{cc}
-w_1 & 0 \\
0 & v_2\end{array}\right) \rho(t) + \left(\begin{array}{c}
w_1\bar{\rho}_1 - F_2 \\
0\end{array}\right) \\
y_V(t) = \left(\begin{array}{c}
F_2 / \rho_1(t)\\
v_2\end{array}\right)\end{array}$\\

\hline
\emph{CC1} &
$\begin{array}{lll}
w_2(\bar{\rho}_2 - \rho_2(t)) & < & v_1\rho_1(t)\\
w_1(\bar{\rho}_1 - \rho_1(t)) & \geq & r_0\\
\rho_2(t) & > & \rho_2^c\end{array}$ &
$\begin{array}{l}
\dot{\rho}(t) = \left(\begin{array}{cc}
0 & w_2 \\
0 & -w_2 \end{array}\right) \rho(t) + \left(\begin{array}{c}
r_0 - w_2\bar{\rho}_2 \\
w_2\bar{\rho}_2 + r_2 - F_2 \end{array}\right) \\
y_f(t) = \left(\begin{array}{cc}
0 & -w_2\\
0 & 0\end{array}\right) \rho(t) + \left(\begin{array}{c}
0\\
F_2\end{array}\right) \\
y_V(t) = \left(\begin{array}{c}
w_2(\bar{\rho}_2 - \rho_2(t)) / \rho_1(t)\\
F_2 / \rho_2(t)\end{array}\right)\end{array}$\\

\hline
\emph{CC2} &
$\begin{array}{lll}
w_1(\bar{\rho}_1 - \rho_1(t)) & < & r_0\\
\rho_2(t) & > & \rho_2^c\end{array}$ &
$\begin{array}{l}
\dot{\rho}(t) = \left(\begin{array}{cc}
-w_1 & w_2 \\
0 & -w_2 \end{array}\right) \rho(t) + \left(\begin{array}{c}
0 \\
w_2\bar{\rho}_2 + r_2 - F_2 \end{array}\right) \\
y_f(t) = \left(\begin{array}{cc}
-w_1 & 0\\
0 & -w_2\end{array}\right)\rho(t) + \left(\begin{array}{c}
w_1\bar{\rho}_1\\
w_2\bar{\rho}_2\end{array}\right) \\
y_V(t) = \left(\begin{array}{c}
w_2(\bar{\rho}_2 - \rho_2(t)) / \rho_1(t)\\
F_2 / \rho_2(t)\end{array}\right)\end{array}$\\

\hline
\end{tabular}
\caption{System (\ref{ctctm0})-(\ref{ctctm3}), (\ref{ctctmoutput0})
for the 2-cell example.}
\label{tab_2cells}
\end{center}\end{table}

\begin{table}[ht]\begin{center}
\begin{tabular}{l|l|l}
Mode & Flow measurement $y_f(t)$ & Velocity measurement $y_V(t)$ \\

\hline
\emph{FF} &
$\begin{array}{lll}
\rho_1(t) & - & \mbox{observable}\\
\rho_2(t) & - & \mbox{observable}\\
\end{array}$ & 
$\begin{array}{lll}
\rho_1(t) & - & \mbox{unobservable}\\
\rho_2(t) & - & \mbox{unobservable}\\
\end{array}$ \\ 

\hline
\emph{FC} &
$\begin{array}{lll}
\rho_1(t) & - & \mbox{observable}\\
\rho_2(t) & - & \mbox{unobservable}\\
\end{array}$ & 
$\begin{array}{lll}
\rho_1(t) & - & \mbox{unobservable}\\
\rho_2(t) & - & \mbox{observable}\\
\end{array}$ \\ 

\hline
\emph{CF1} &
$\begin{array}{lll}
\rho_1(t) & - & \mbox{unobservable}\\
\rho_2(t) & - & \mbox{observable}\\
\end{array}$ & 
$\begin{array}{lll}
\rho_1(t) & - & \mbox{observable}\\
\rho_2(t) & - & \mbox{unobservable}\\
\end{array}$ \\ 

\hline
\emph{CF2} &
$\begin{array}{lll}
\rho_1(t) & - & \mbox{observable}\\
\rho_2(t) & - & \mbox{observable}\\
\end{array}$ & 
$\begin{array}{lll}
\rho_1(t) & - & \mbox{observable}\\
\rho_2(t) & - & \mbox{unobservable}\\
\end{array}$ \\ 

\hline
\emph{CC1} &
$\begin{array}{lll}
\rho_1(t) & - & \mbox{unobservable}\\
\rho_2(t) & - & \mbox{observable}\\
\end{array}$ & 
$\begin{array}{lll}
\rho_1(t) & - & \mbox{observable}\\
\rho_2(t) & - & \mbox{observable}\\
\end{array}$ \\ 

\hline
\emph{CC2} &
$\begin{array}{lll}
\rho_1(t) & - & \mbox{observable}\\
\rho_2(t) & - & \mbox{observable}\\
\end{array}$ & 
$\begin{array}{lll}
\rho_1(t) & - & \mbox{observable}\\
\rho_2(t) & - & \mbox{observable}\\
\end{array}$ \\ 

\hline
\end{tabular}
\caption{Observability properties of the 2-cell system.}
\label{tab_2cellobservability}
\end{center}\end{table}

\begin{figure}[htbp]
\centerline{
\includegraphics[width=4in]{figures/pp2cells.eps}}
\caption{Phase portrait of the 2-cell system and the partition into regions
with affine dynamics.
Green and red dots together with yellow and orange line segments depict
the set of equilibria for this system.}
\label{figpp2cells}
\end{figure}

Figure \ref{figpp2cells} illustrates the dynamics of the 2-cell system and
shows how the state space is partitioned into six regions \emph{FF}, \emph{FC},
\emph{CF1}, \emph{CF2}, \emph{CC1} and \emph{CC2}, in each of which the
system dynamics is affine and the observation equation for $y_f(t)$ is
affine\footnote{The number of affine modes grows exponentially with the
dimension of the system.
For general $N$ the number of affine modes is $2^N+2$.}.
From the observation equations in Table \ref{tab_2cells} it is easy to see
which of part of the state $\rho(t)$ can be determined from which measurements
in each affine mode.

Table \ref{tab_2cellobservability} summarizes the observability
properties of the 2-cell system.
We can see that measuring the flow $y_f(t)$ is sufficient for estimating
the density only in the modes \emph{FF}, \emph{CF2} and \emph{CC2},
whereas measuring velocity $y_V(t)$ is sufficient for estimating
the density only in the modes \emph{CC1} and \emph{CC2}.
In modes \emph{FC} and \emph{CF1} it is not enough to mesure just flow
or just velocity to determine density $\rho(t)$ --- here we need both
measurements, $y_f(t)$ and $y_V(t)$.

In \cite{munoz03} the authors investigate the observability of CTM
for the case of flow measurement and argue that using piecewise affine
approximation of the original system with just two affine modes, which
correspond to our \emph{FF} and \emph{CC2} switching randomly between those,
is fair enough for the purpose of density estimation.
The obvious setback of such approach is that this approximation of CTM is
no longer a conservation law and thus has different dynamics.
The reason why estimation with the two-mode switching model produces decent
results is in the fact that the demand is practically never feasible --- it is
either strictly feasible or infeasible.
It was shown in \cite{gomes08} that for both, strictly feasible and
infeasible demand, the system has globally asymptotically stable equlibrium
points, and thus, it has to end up in either all free flow \emph{FF}
or all congested \emph{CC2} mode.










\subsection{Problem Setting}\label{subsec_problemsetting}
In order to find bottlenecks and based on their locations partition the
freeway into segments, each of which could be treated separately using
the valuable monotonicity property of CTM, we need to have constant
off-ramp split ratios and constant demand.
Off-ramp split ratios $\beta_i$, $i=1..N$, are assumed to be known.
We shall also assume the demand to be known with some uncertainty
and piecewise constant: $r\sim\NN(\bar{r}, \Sigma_r)$, where the mean
value lies in the box defined by $r_i^-\leq r_i\leq r_i^+$, $i=1..N$
($r_i^-$ and $r_i^+$ are known), and $\Sigma_r$ is known covariance matrix.
It can be updated every given time period $T$
(for example, $T=1/4$ hour)\footnote{The problem of demand prediction
can be addressed by various available approaches,
such as \cite{lin01, okutani84, smith02}.
Systems such as PeMS \cite{pems} provide historic measurement data.}.

Cell capacities $F_i$ are also known with uncertainty.
Established in the process of calibration they are likely to be estimated
with error.
Therefore, we shall assume $F\sim\NN(\bar{F}, \Sigma_F)$, where the mean value
$\bar{F}$ lies in the box defined by $F_i^-\leq \bar{F}_i\leq F_i^+$, $i=1..N$
($F_i^-$ and $F_i^+$ are known), and $\Sigma_F$ is known covariance matrix.
Accordingly, we define
\begin{equation}
\rho_i^{c-} = \frac{F_i^-}{v_i} \; \mbox{ and } \;
\bar{\rho}_i^- = \frac{F_i^-}{v_i} + \frac{F_i^-}{w_i},
\label{rhobarminus}
\end{equation}
and
\begin{equation}
\rho_i^{c+} = \frac{F_i^+}{v_i} \; \mbox{ and } \;
\bar{\rho}_i^+ = \frac{F_i^+}{v_i} + \frac{F_i^+}{w_i}.
\label{rhobarplus}
\end{equation}

Define $f_i^+$, $i=0..N$, as in (\ref{2d})-(\ref{2a}) with $r_i$ substituted
by $r_i^+$.
Then, cell $i$ is considered to be a bottleneck if the equilibrium flow $f_i$
satisfies $f_i^+\geq F_i^-$.

Denoting the right hand side of system (\ref{ctctm0}) as $g(t, \rho(t))$ and
adding parametric uncertainty $\xi(t)$ that comes from demand and capacity,
we substitute this system with
\begin{equation}
\dot{\rho}(t) = g(t, \rho(t)) + \xi(t).\label{ctctmstate}
\end{equation}
Here we assume $\xi(t)\sim\NN(\bar{\xi}, \Sigma_\xi)$, where the mean
value $\bar{\xi}$ lies in the box
$-\xi_i^0\leq\bar{\xi}_i\leq\xi_i^0$, $i=1..N$,
with known $\xi_i^0\geq 0$, and $\Sigma_\xi$ is known covariance matrix.

Some measurement noise $\omega(t)$ is present.
We assume $\omega(t) \sim \NN(\bar{\omega}, \Sigma_\omega)$,
where the mean value $\bar{\omega}$ lies in the box
$-\omega_i^0\leq \bar{\omega}_i\leq \omega_i^0$ for every $i$,
with known $\omega_i^0\geq 0$, and $\Sigma_\omega$ is known covariance matrix.
So, the observation equation (\ref{ctctmoutput0}) is replaced with
\begin{equation}
y(t) = h(t, \rho(t)) + \omega(t).\label{ctctmoutput}
\end{equation}
To ensure complete observability of the system, we shall assume that both,
flow measurment $y_f(t)$ and velocity measurement $y_V(t)$ are
available\footnote{Had there be no noise ($\omega(t)=0$),
we could determine density $\rho(t)$
without using the state equation (\ref{ctctmstate}), only from the
measurements $y_f(t)$ and $y_V(t)$ and formula (\ref{density-flow-speed}):
\[ \rho_i(t) = \frac{y_{fi}(t)}{y_{Vi}(t)}, \;\;\;
i=1..N .\label{rho_measured} \]
In reality, however, the noise is present.}.
We shall describe the state estimation for the system
(\ref{ctctmstate}), (\ref{ctctmoutput}) with
$y(t) = \left(\begin{array}{c}
y_f(t)\\
y_V(t)\end{array}\right)$.





\newcommand{\hrho}{\hat{\rho}}
\newcommand{\hy}{\hat{y}}




\subsection{Set-Valued Estimation}\label{subsec_setvalued}
First, we consider the non-stochastic case.
That is, covariance matrices $\Sigma_\xi = \Sigma_\omega \equiv 0$.
So, the uncertainty and the measurement noise belong to box-sets that may
evolve in time:
\begin{equation}
\xi(t) \in \left\{\bar{\xi} \in {\bf R}^N ~|~
-\xi_i^0(t) \leq \bar{\xi}_i \leq \xi_i^0(t), \; i=1..N \right\},
\label{xi_bounds}
\end{equation}
where $\xi^0(t) = \left(\begin{array}{c}
\xi_1^0(t)\\
\vdots \\
\xi_N^0(t)\end{array}\right)$, $0\leq t\leq T$, is known nonnegative function;
and 
\begin{eqnarray}
\omega_f(t) & \in & \left\{\bar{\omega}_f \in {\bf R}^N ~|~
-\omega_{fi}^0(t) \leq \bar{\omega}_{fi} \leq \omega_{fi}^0(t),
\; i=1..N \right\}, \label{omega_bounds_f} \\
\omega_V(t) & \in & \left\{\bar{\omega}_V \in {\bf R}^N ~|~
-\omega_{Vi}^0(t) \leq \bar{\omega}_{Vi} \leq \omega_{Vi}^0(t),
\; i=1..N \right\}, \label{omega_bounds_V}
\end{eqnarray}
where subscript $_f$ indicates values related to flow measurements,
subscript $_V$ indicates values related to velocity measurements,
and $\omega_f^0(t) = \left(\begin{array}{c}
\omega_{f1}^0(t)\\
\vdots \\
\omega_{fN}^0(t)\end{array}\right)$, $\omega_V^0(t) = \left(\begin{array}{c}
\omega_{V1}^0(t)\\
\vdots \\
\omega_{VN}^0(t)\end{array}\right)$, $0\leq t\leq T$,
are known nonnegative functions.

Denote (\ref{ctctmstate}$^-$) system (\ref{ctctmstate}) with
$\xi(\cdot) = -\xi^0(\cdot)$, and (\ref{ctctmstate}$^+$) --- system
(\ref{ctctmstate}) with $\xi(\cdot)=\xi^0(\cdot)$.

Data measurements are collected at times
$0 = \tau_0 < \tau_1 < \cdots < \tau_{k-1} < \tau_k = T$.
The state estimation algorithm works is described below.

\begin{enumerate}
\item First, we must establish the initial conditions --- that is, determine
bounds for the initial state $\rho(0)$.
At time $t=0$ we take the first measurement
$\hy(0) = \left(\begin{array}{c}
\hy_f(0)\\
\hy_V(0)\end{array}\right)$, which together with (\ref{ctctmoutput}) gives us
\[ \hrho(0) = \left(\begin{array}{c}
\frac{\hy_{f1}(0) - \omega_{f1}(0)}{\bar{\beta}_1(\hy_{V1}(0) - \omega_{V1}(0))} \\
\vdots \\
\frac{\hy_{fN}(0) - \omega_{fN}(0)}{\bar{\beta}_N(\hy_{VN}(0) - \omega_{VN}(0))}
\end{array}\right), \]
where $\omega_{fi}(0)$, $\omega_{Vi}(0)$, $i=1..N$,
are known each within interval
$\omega_{fi}(0)\in[-\omega_{fi}^0(0), \omega_{fi}^0(0)]$ and
$\omega_{Vi}(0)\in[-\omega_{Vi}^0(0), \omega_{Vi}^0(0)]$.
So,
\begin{equation}
\hrho_i(0) \in [\hrho_i^-(0), \hrho_i^+(0)], \;\;\; i=1..N. \label{hrho0}
\end{equation}

Suppose, the following conditions hold:
\begin{equation}
\hy_{fi}(0) - \omega_{fi}^0(0) > 0,
\label{cond1}
\end{equation}
and 
\begin{equation}
\hy_{Vi}(0) - \omega_{Vi}^0(0) > 0.
\label{cond2}
\end{equation}
Then, following (\ref{int_plus}) and (\ref{int_quotient1})
from Appendix \ref{app_intervals}, we get
\begin{equation}
\hrho_i^-(0) =
\frac{\hy_{fi}(0) - \omega_{fi}^0(0)}{\bar{\beta}_i(\hy_{Vi}(0) + \omega_{Vi}^0(0))} \;
\mbox{ and } \;
\hrho_i^+(0) =
\frac{\hy_{fi}(0) + \omega_{fi}^0(0)}{\bar{\beta}_i(\hy_{Vi}(0) - \omega_{Vi}^0(0))}.
\label{hinterval0}
\end{equation}
If condition (\ref{cond2}) does not hold, set
\begin{equation}
\hrho_i^-(0) = \hrho_i^+(0) = \bar{\rho}_i^+, \label{hinterval0a}
\end{equation}
where $\bar{\rho}_i^+$ is defined in (\ref{rhobarplus}),
and if condition (\ref{cond2}) holds but condition (\ref{cond1}) does not, set
\begin{equation}
\hrho_i^-(0) = \hrho_i^+(0) = 0. \label{hinterval0b}
\end{equation}
We may have some apriori knowledge about the bounds of the initial
state\footnote{Information about the initial state may come from the
previous estimation.}:
\begin{equation}
\rho_i(0) \in [\rho_i^-(0), \rho_i^+(0)], \;\;\; i=1..N, \label{rho0}
\end{equation}
then, the initial state should belong to the intersection of sets
(\ref{hrho0}) and (\ref{rho0}):
\begin{equation}
\left\{\rho \in {\bf R}_+^N ~|~
\rho_i \in [\rho_i^-(0), \rho_i^+(0)] \cap [\hrho_i^-(0), \hrho_i^+(0)] ~
i=1..N\right\} . \label{initialset}
\end{equation}
If this set is not empty, update the initial state bounds:
\begin{equation}
\rho_i^-(0) = \max(\rho_i^-(0), \hrho_i^-(0)) \; \mbox{ and } \;
\rho_i^+(0) = \min(\rho_i^+(0), \hrho_i^+(0)) \;\;\; i=1..N,
\label{rhobounds0a}
\end{equation}
else, if (\ref{initialset}) is empty,
or there is no apriori information about the initial state, set
\begin{equation}
\rho_i^-(0) = \hrho_i^-(0) \; \mbox{ and } \;
\rho_i^+(0) = \hrho_i^+(0) \;\;\; i=1..N,
\label{rhobounds0b}
\end{equation}
Finally, we make sure that our density bounds do not exceed jam density
by setting
\begin{equation}
\rho_i^-(0) = \min(\rho_i^-(0), \bar{\rho}_i^+) \; \mbox{ and } \;
\rho_i^+(0) = \min(\rho_i^+(0), \bar{\rho}_i^+) \;\;\; i=1..N,
\label{rhobounds0c}
\end{equation}
Obviously,
\[ \rho^-(0) = \left(\begin{array}{c}
\rho_1^-(0)\\
\vdots \\
\rho_N^-(0)\end{array}\right) \leq \left(\begin{array}{c}
\rho_1^+(0)\\
\vdots \\
\rho_N^+(0)\end{array}\right) = \rho^+(0) .\]



Now, for $j = 1..k$ repeat the following steps.

\item Denote $\rho^-(\tau)$ and $\rho^+(\tau)$,
$\tau_{j-1}\leq\tau\leq\tau_{j}$,
the trajectories of systems (\ref{ctctmstate}$^-$) and (\ref{ctctmstate}$^+$)
with initial states $\rho^-(\tau_{j-1})$ and $\rho^+(\tau_{j-1})$ accordingly.
Due to the fact that system (\ref{ctctmstate}) is monotonic,
it is true that
\[ \rho^-(\tau_{j-1})\leq\rho^+(\tau_{j-1}) \Rightarrow
\rho^-(\tau)\leq\rho^+(\tau) \]
for all $\tau\geq\tau_{j-1}$.
Hence, any trajectory $\rho(\tau)$, $\tau_{j-1}\leq\tau\leq\tau_j$,
of system (\ref{ctctmstate}) with initial state $\rho(\tau_{j-1})$ such that
$\rho^-(\tau_{j-1})\leq\rho(\tau_{j-1})\leq\rho^+(\tau_{j-1})$ satisfies
\begin{equation}
\rho^-(\tau)\leq\rho(\tau)\leq\rho^+(\tau), \;\;\; \tau_{j-1}\leq\tau\leq\tau_j.
\label{monotonicity}
\end{equation}
Thus, after computing trajectories $\rho^-(\tau)$ and $\rho^+(\tau)$ for
$\tau_{j-1}\leq\tau\leq\tau_j$, we can state that for any $\tau$ in this
interval, state $\rho(\tau)$ of system (\ref{ctctmstate})
satisfies the inclusion 
\begin{equation}
\rho(\tau) \in \left\{\rho\in{\bf R}^N_+ ~|~
\min(\rho_i^-(\tau),\bar{\rho}_i^+)
\leq\rho_i(\tau)\leq\min(\rho_i^+(\tau), \bar{\rho}_i^+), ~
i=1..N\right\}
\label{monotonicity2}
\end{equation}

\item At time $t=\tau_j$ we receive the measurement data
$\hy(\tau_j) = \left(\begin{array}{c}
\hy_f(\tau_j)\\
\hy_V(\tau_j)\end{array}\right)$, which, like in step 1, gives us the
estimate
\begin{equation}
\hrho(\tau_j) \in \left\{\hrho\in{\bf R}^N_+ ~|~
\hrho_i^-(\tau_j)\leq\hrho(\tau_j)\leq\hrho_i^+(\tau_j), ~ i=1..N \right\},
\label{rhoestimate}
\end{equation}
where
\begin{equation}
\hrho_i^-(\tau_j) =
\frac{\hy_{fi}(\tau_j) - \omega_{fi}^0(\tau_j)}{\bar{\beta}_i(\hy_{Vi}(\tau_j)
+ \omega_{Vi}^0(\tau_j))} \;
\mbox{ and } \;
\hrho_i^+(\tau_j) =
\frac{\hy_{fi}(\tau_j) + \omega_{fi}^0(\tau_j)}{\bar{\beta}_i(\hy_{Vi}(\tau_j)
- \omega_{Vi}^0(\tau_j))}
\label{hinterval}
\end{equation}
if conditions
\begin{equation}
\hy_{fi}(\tau_j) - \omega_{fi}^0(\tau_j) > 0
\label{cond11}
\end{equation}
and 
\begin{equation}
\hy_{Vi}(\tau_j) - \omega_{Vi}^0(\tau_j) > 0
\label{cond22}
\end{equation}
are satisfied.
Otherwise, if condition (\ref{cond22}) does not hold, set
\begin{equation}
\hrho_i^-(\tau_j) = \hrho_i^+(\tau_j) = \bar{\rho}_i^+, \label{hintervala}
\end{equation}
or, if (\ref{cond22}) holds, but (\ref{cond11}) does not, set
\begin{equation}
\hrho_i^-(\tau_j) = \hrho_i^+(\tau_j) = 0. \label{hintervalb}
\end{equation}
Now we perform the correction in our estimation by taking intersection
of the box-sets (\ref{monotonicity2}) and (\ref{rhoestimate}):
\begin{eqnarray*}
&& \left\{\rho\in{\bf R}^N_+ ~|~
\min(\rho_i^-(\tau_j),\bar{\rho}_i^+)
\leq\rho_i(\tau_j)\leq\min(\rho_i^+(\tau_j), \bar{\rho_j}_i^+), ~
i=1..N \right\} \\
&& \cap \left\{\hrho\in{\bf R}^N_+ ~|~
\hrho_i^-(\tau_j)\leq\hrho(\tau_j)\leq\hrho_i^+(\tau_j), ~ i=1..N \right\}.
\end{eqnarray*}
If this intersection is empty, set
\begin{equation}
\rho^-(\tau_j) = \hrho^-(\tau_j) \; \mbox{ and } \;
\rho^+(\tau_j) = \hrho^+(\tau_j).
\label{boundsupdate2}
\end{equation}
Update the bounds $\rho^-(\tau_j)$ and
$\rho^+(\tau_j)$:
\begin{equation}
\rho^-(\tau_j) = \left(\begin{array}{c}
\max(0, \hrho_1^-(\tau_j), \min(\rho_1^-(\tau_j), \bar{\rho}_1^+))\\
\vdots\\
\max(0, \hrho_N^-(\tau_j), \min(\rho_N^-(\tau_j), \bar{\rho}_N^+))
\end{array}\right)
\label{boundsupdate1a}
\end{equation}
and
\begin{equation}
\rho^+(\tau_j) = \left(\begin{array}{c}
\max(0, \min(\hrho_1^+(\tau_j), \rho_1^+(\tau_j), \bar{\rho}_1^+))\\
\vdots\\
\max(0, \min(\hrho_N^+(\tau_j), \rho_N^+(\tau_j), \bar{\rho}_N^+))
\end{array}\right).
\label{boundsupdate1b}
\end{equation}

\item Set $j=j+1$, and while $j\leq k$ repeat steps 2, 3 and 4.
\end{enumerate}









